Algorithm for Linear Response Functions at Finite Temperatures: 
Application to ESR spectrum of s = | Antiferromagnet Cu benzoate 
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We introduce an efficient and numerically stable method for calculating linear response functions 
x(q, of quantum systems at finite temperatures. The method is a combination of numerical 
solution of the time-dependent Schroedinger equation, random vector representation of trace, and 
Chebyshev polynomial expansion of Boltzmann operator. This method should be very useful for a 
wide range of strongly correlated quantum systems at finite temperatures. We present an application 
to the ESR spectrum of s = | antiferromagnet Cu benzoate. 
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PACS numbers: 75.10.Jm, 75.40.Mg 



I. INTRODUCTION 

Computational physicists often face to problems of cal- 
culating the linear response functions or the real-time 
Green's functions of quantum manybody systems with 
degree of freedom N ~ 10 6 or more. Direct diagonal- 
ization of the Hamiltonian matrix is the most naive and 
powerful method for modest system size N < 10 3 . How- 
ever, it becomes prohibitive for larger system because the 
computational time grows as N 3 . Therefore efficient nu- 
merical algorithms, such as Quantum Monte Carlo meth- 
ods, Lanczos Methods, and Kernel Polynomial Methods 
have been developed and applied to various problems. 

Quantum Monte Carlo methods [jl], |) can generate the 
Green's functions of very large systems since QMC does 
not need to store the wave functions. They have been suc- 
cessfully used for evaluating the imaginary-time Green's 
functions and, therefore, various thermodynamic quanti- 
ties. However, for evaluating dynamical quantities such 
as AC conductivity, one has to rely on numerical analytic 
continuation (e.g. Maximum Entropy Method [p| ) to ob- 
tain the real-time Green's function from the imaginary- 
time one. This procedure is not straightforward because 
the statistical errors are amplified by numerical analyt- 
ical continuation and the default model for the MEM 
must be assumed a priori. 

Lanczos Methods [Q, || have been useful techniques for 
evaluating dynamical responses of relatively large sys- 
tems. The LM uses Lanczos recursion formula with Ma- 
trix Vector Multiplications (\<p') = H\(f>)) to tridiagonal- 
ize the Hamiltonian matrix, which leads to a continued 
fraction representation of the real-time Green's function. 
The drawback of LM is numerical instability of Lanczos 
recursion formula caused by accumulation of roundoff er- 
rors for large numbers of MVM's. Recently, LM was ex- 
tended to finite temperatures (Finite Temperature LM) 
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by introducing random sampling over the ground and ex- 
cited eigenstates ||. However, FTLM has a weak point 
that the number of excited eigenstates to be calculated 
increases rapidly as temperature or system size becomes 
large. Therefore reduction of computational costs by ex- 
ploiting symmetries of the system becomes crucial. 

Kernel Polynomial Method [0, || calculates density of 
states and linear response functions by using Chebyshev 
polynomial expansion ||, |l0|, [ll]] of broadened delta func- 
tions. The Chebyshev polynomials are obtained through 
Chebyshev recursion formula by using MVM's. Unlike 
Lanczos recursion, Chebyshev recursion avoids accumu- 
lation of numerical roundoff errors even for large numbers 
of MVM's. The use of KPM for linear response functions 
has been limited to the ground state calculation. 

Fermi- Weighted Time-Dependent Method jl2| is a 
combination of conventional time-dependent method [13 , 
O, |li|, random vector representation of trace fllq , 
and Chebyshev polynomial expansion of step function 
0(Ef — H) to extract occupied and unoccupied one- 
particle states of noninteracting fermi particles in the 
ground state. This method was successfully applied to 
optical absorption of silicon nano-crystallitesflj. 

Our new algorithm, the Boltzmann- Weighted Time- 
Dependent Method, calculates linear response functions 
at finite temperatures by using the Boltzmann-weight 
operator in place of the Fermi-weight operator. The 
BWTDM has an advantage to FTLM that the compu- 
tational time does not increase as temperature increases, 
because BWTDM does not calculate any eigenstates at 
all. Therefore BWTDM does not need use of symmetries 
and can easily treat disordered systems. 

In this Letter, we introduce BWTDM and ap- 
ply it to study Electron Spin Resonance spectrum 
of one-dimensional s = \ antiferromagnet Cu ben- 
zoate Cu(CqH 5 COO)2 ■ 3H 2 0, especially the dynamical 
crossover between spinon excitation and breather exci- 
tation as a function of temperature. Nonperturbative 
analytical calculation of this phenomena is a very diffi- 
cult problem which still open. Even within Sine Gor- 
don theory, the linear response functions have been stud- 
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ied only at high temperature regime (by perturbation) 
and at zero temperature. We compare our numerical re- 
sults with preceding experimental and theoretical studies 
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II. THE MODEL 

The effective Hamiltonian of one-dimensional s — 1 
antiferromagnet Cu benzoate may be written p5| , ^61 |27 
as 



H = 



E 



[Jsj ■ s j+ i - g^ B H x s* + (-l) j giJ, B hsV] (.1) 



The first term is the isotropic Heisenberg antiferromag- 
net where Sj is spin operator at the j-th site, J/k B 
1 7.2 1\ = lA8(meV) is the exchange interaction Jl8| , 
and the sums are taken over N sp i n spins with periodic 
boundary conditions. The second is the normal Zee- 
man term due to the external uniform magnetic field 
H x . The third is the staggered Zeeman term due to 
induced staggered magnetic field h = 0.065H X originat- 
ing from the Dzyaloshinskii-Moriya (DM) interaction and 
the staggered component of the g tensor ||]|. In the fol- 
lowing we assume that g = 2.25(2.29) with H x along 
c(c")-direction. Unit of magnetic field and frequency 
become J I g c \i B = H-4(T), J/g c »HB = H-2(T), and 
J/h = 359(GHz), respectively. 

According to the linear response theory, absorption of 
electromagnetic waves of frequency uj and polarization (i 
is expressed as 



(2) 



where Hr is the amplitude of the electromagnetic waves 
and Xuu(9) w ) i s the imaginary part of the dynamical 
magnetic susceptibility, 

/>OC 

x'^tJ) = (1 - e-P») Im lim / dte'^-^g^t) 

(3) 

at inverse temperature (3 = l/(k B T). The correlation 
function is defined by 

g${t) = Tr [e-^ H M^ q e +im Ml q e- iHt ] /Tr [e~ pB ] (4) 

where the magnetization operator is = 

^' =1 S ^"V^' 



III. ALGORITHM 

The essence of BWTDM is evaluation of eq. (|j) by 
using Boltzmann-weighted random vectors, \$ Bo itz) — 
e"^/ 2 |$), 

, (f] (( (*Boit*\M!t q e+ iHt \M!t q \e- iHt \* Boltz ) )) 

9q 1 J (( ($Boltz\*Boltz) » ^ j 
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FIG. 1: ESR spectra of normal polarization, I z (lo), calculated 
with P = 2 ~ 16, H x = 1.0, iV 3 pi„ = 16, n = 0.01 and 
N ran d = 16. S and Bi stand for spinon excitation and first 
breather excitation, respectively. 



In the above equation, the uniform random vector, |$), 
is defined by 



i*) = E \ n )^ 



(6) 



where {\n}} is the basis set composed of eigenvectors of 
Sj's, and are a set of complex random variables satis- 
fying 



(7) 



Here (( • )) stands for the statistical average. Then 
(( )) gives the trace of X, 

((($|X|$>» =J2((CCn))(n\X\n) 

n 

+ E« &z» » 

= 1 £{n\X\n)=tr[X\ (8) 

n 

where we used eq. (^) . The Boltzmann operator is evalu- 
ated by the Chebyshev polynomial expansion || [lo[ [ll] . 



e- TH \d>) =^k{r)T k {H)\4>) 



(9) 



where afc(r) is the coefficient of the Chebyshev poly- 
nomial expansion, and \<j>) is an arbitrary vector. The 
Chebyshev polynomial of k-th order T k (x) is calculated 
by the Chebyshev recursion formula, 

T k+1 {H)\4>) = 2HT k {H)\4>) - T).-i(H)\<fi). (10) 

The time-evolution operator \<j>;t + At) — e~ lHAt \4>;t} 
can be evaluated with numerical stability by the leap 
frog method Q , 

\<f>;t + At) = -2iHAt\<j>;t) + \<f>;t- A). (11) 

In real calculation, we first generate N complex ran- 
dom numbers £ n and construct |<fe) according to eq. ra). 
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FIG. 2: (a) Peak width (FWHM) of spinon (S) and first 
breather (Si) excitations calculated with /3 = 0.5 ~ 9.5, H x — 
1.0, N S pi n = 16, r) = 0.01 and N ran d = 16. The horizontal 
solid line represents the width due to the finite resolution 
rj = 0.01. 



The \$Boltz) can be computed with numerical stability by 
applying eq. @ repeatedly to |$). Then ($ Bo it *|$ 'Bolt z) 
gives a sample of the denominator of eq. (||). Here, 
we introduce another vector \$m^ ) — M+ q \&Boitz) and 
calculate the time evolution of \^Boitz) and \&m^ ) by 
cq.([ll]). At each time t = nAt, the matrix element 
($ m m ; t\M+ q \<& Boitz] t) gives a sample of numerator of 
eq. (^). After calculating the denominator and numera- 
tor for N rani i random vectors, the ratio of their averages 
gives g%(t) of eq. (||). Then Xuuili w ) with a finite fre- 
quency resolution r\ is calculated by using Gaussian filter, 



X " (q,u) = (1 - e-**) Im / dte**£{t) 



Ae-^ t2 / 2 



(12) 

where T max is related to 77 by T ma:E ~ 77 1 . In order 
to avoid finite size effect of spinon excitations, T max is 
chosen so that T max < N spin /v spin where v spin ~ 1 is 
the spin wave velocity |27|| . This limits the finest fre- 
quency resolution. Note that this is not the limit due 
to the algorithm but due to the finite size model. Much 
finer resolution can be used for breather modes, which 
are spatially localized. 

Two more numerical details: First, the fluctuation in 
the result is inversely proportional to the square root of 
the number of eigenstates participating in \^Boitz) mul- 
tiplied by N ranc i [0 • Therefore systems at lower temper- 
atures require more random vectors for the same quality 
of the result. Second, since the Chcbyshev polynomials 
are defined only in the range [—1,-1-1], the Hamiltonian 
should be rescaled when eq. (||) is applied, so that all 
eigenvalues of H lie in the range [—1, +1]. 



IV. RESULTS AND DISCUSSIONS 

Figure 1 shows the ESR spectra of normal polarization 
(Faraday configuration) I z (lu) calculated with various 
temperatures, and Fig. 2 shows the width and peak fre- 
quency of the spinon excitation S and the first breather 
excitation B\ obtained by fitting two Gaussian peaks to 



FIG. 2: (b) Peak frequency of spinon (5*) and first breather 
(B\) excitations calculated with j3 = 0.5 ~ 9.5, H x — 1.0, 
N spin = 16, 7/ = 0.01 and N ran d = 16. 



the calculated spectra with the least square method. At 
high temperatures [ft < 5), S is outstanding in the spec- 
trum, and its peak shifts to higher frequency and becomes 
broadened as temperature decreases while B\ becomes 
narrower. At low temperatures (/3 > 5), B\ prevails 
while its peak frequency is almost constant and its width 
becomes much smaller than the numerical resolution rj. 
This crossover behavior between spinon excitation and 
first breather excitation has become computable by the 
invention of BWTDM, and the result is consistent with 
both experimental Jl9| and field-theoretical results [p3| . 
In addition, we can observe in Fig. 2 the shift of B\ at 
high temperatures and the narrowing of S at low tem- 
peratures. We believe these tendencies are real physics 
but not numerical artifact though further experimental 
and theoretical studies are necessary to confirm it. Sev- 
eral other weak peaks appear in our results as well as in 
experimental results, which are supposed to be higher- 
order breathers and transitions between excited states. 
However, we are not going to analyze them here. 

Figure 3 shows the induced energy gap E (H X ) calcu- 
lated from the Bi frequency u> by solving pll 



frw = J(g[i B H x ) 2 + {E g (H x )y 



For H < 1, E g (H x ) scales to H x ^ 3 as predicted by Os- 
hikawa and Affleck |25). For 1 < H x < 4, our result seems 

2 /3 

to deviate from the H x -scaling and have a dip around 
H x ' ~ 2. This seems consistent with the experimen- 
tal results and the density-matrix renormalization-group 
study |Q . At such a high field, assumption of Sine Gor- 
don theory, that the H x is weak and the frequency is 
small, is probably broken. Since the gap is calculated as 
a small difference of two large numbers, fkv and gj.isH x , 
the agreement of our result with the experimental value 

2/3 

and the H x -scaling at low fields demonstrates the nu- 
merical accuracy of our algorithm. 

In summary, we developed an efficient and stable al- 
gorithm for linear response functions at finite temper- 
atures, and studied the ESR spectra of s — \ antifcr- 
romagnet Cu benzoate for a wide range of temperature 
and magnetic field. We reproduced experimental results 
of the spinon-breather crossover as a function of temper- 
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FIG. 3: Solid curve with filled squares shows E g (H x ) calcu- 
lated from our results with j3 = 16, N Bp i n = 16, r\ = 0.01 and 
N ran d = 16. Solid curve with filled circles shows E g (H x ) for 
H x || c — axis calculated from the experimental B\ peak read 
from Figure 4 of Ref. . Thin line shows a&tcj — cH^ 3 to 
the experimental result. [ p^[ . 

ature. Temperature dependence of the width and shift of 
the peaks are also calculated, which are consistent with 
experiments and analytical theories. The calculated fre- 
quency of B\ as a function of H x agrees well with both ex- 
perimental and field-theoretical results at low fields, and 
reproduces the deviation of experimental results from the 
Sine Gordon theory at high fields, where the low energy 



assumption of the Sine Gordon theory may be broken 
and the choice of the compactification radius R becomes 
ambiguous ]24|. The advantage of BWTDM is being ap- 
plicable to finite temperatures, strong magnetic field and 
high frequency while its weak point is finite size effects 
(N spin < 20). The computational cost of BWTDM is 
moderate. Calculation of one curve in Fig. 1, for ex- 
ample, requires approximately 30 minutes with 8 CPU's 
of Fujitsu VPP5000 vector-parallel computer. We hope 
that this algorithm stimulates further numerical inves- 
tigation of dynamical properties of quantum manybody 
systems at finite temperatures in general. 
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